AquaMaps ingest

Initialize

• AquaMaps Environmental Dataset: Half-Degree Cells Authority File (HCAF)

library(tidyverse)
library(glue)
library(fs)
library(raster)
library(leaflet)

# dir & csv's
#dir_dat   <- "/Users/bbest/Desktop/aquamaps_ver0816c"
dir_dat   <- "/Volumes/Best HD/mbon_data_big/aquamaps"
spp_csv   <- file.path(dir_dat, "speciesoccursum_ver0816c.csv")
obs_csv   <- file.path(dir_dat, "occurrencecells_ver0816c_head-100K-rows.csv")
#cells_csv <- file.path(dir_dat, "hcaf_species_native_ver0816c_head-100K-rows.csv")
cells_csv <- file.path(dir_dat, "hcaf_species_native_ver0816c.csv")
env_csv   <- file.path(dir_dat, "hcaf_v6.csv")
spp_csv   <- file.path(dir_dat, "spp.csv")
dir_spp   <- file.path(dir_dat, "species_rasters")

spp <- read_csv(spp_csv) # head(occ)
#cells <- read_csv(cells_csv)
#obs   <- read_csv(obs_csv)

# AquaMaps raster specifications for 0.5 degree global raster
r_na <- raster(
  xmn = -180, xmx = 180, ymn = -90, ymx = 90, 
  resolution=0.5, crs=leaflet:::epsg4326)
r_a <- area(r_na)

r_1 <- setValues(r_na, 1)

# writeRaster(r_a, path(dir_dat, "global_area.tif"))
# writeRaster(r_1, path(dir_dat, "global_1.tif"))

env <- read_delim(env_csv, delim = "\t")
# TODO: use env.OceanArea for r_a

raster_trim <- function(r){
  # source: /Users/bbest/github/nrelutils/R/raster_utils.R
  
  # raster::trim() crazy slow compared to this
  r_m          <- is.na(raster::as.matrix(r))
  col_notna    <- which(colSums(r_m) != nrow(r))
  row_notna    <- which(rowSums(r_m) != ncol(r))
  extent_notna <- raster::extent(
    r,
    row_notna[1], row_notna[length(row_notna)],
    col_notna[1], col_notna[length(col_notna)])
  raster::crop(r, extent_notna)
}

sp_id_to_tif <- function(sp_id){
  # , cells=cells, r_na=r_na, dir_spp=dir_spp
  
  # get cells
  cells_sp <- cells %>%
    filter(SpeciesID == sp_id) %>%
    mutate(
      cell = cellFromXY(r_na, matrix(data = c(CenterLong, CenterLat), ncol=2)))
  
  # assign to raster
  r_sp <- r_na
  r_sp[cells_sp$cell] <- cells_sp$probability
  r_sp <- raster_trim(r_sp) # plot(r_sp)
  
  # write out
  r_sp_tif <- glue("{dir_spp}/{sp_id}.tif")
  writeRaster(r_sp, r_sp_tif, overwrite=T)
}

# calculate area of raster cells (km2)
r_a <- area(r_na)

# summarize and plot
summary(values(r_a))
plot(r_a)

Table (*.csv) to Map (*.tif) for All Species

cells <- read_csv(cells_csv)

spp_ids <- distinct(cells, SpeciesID) %>% pull(SpeciesID)

# takes 8.8 hours, 472 MB
for (i in 1:length(spp_ids)){ # i = 1
  cat(sprintf("%04d of %d: %s -- %s\n", i, length(spp_ids), sp_id, Sys.time()))
  # https://www.aquamaps.org/preMap2.php?cache=1&SpecID=Fis-22812
  sp_id <- spp_ids[i]
  
  sp_id_to_tif(sp_id)
}

# estimate time to completion ----
  
# 0001 of 24904: Fis-22812 -- 2018-07-01 20:49:01
# 0349 of 24904: Fis-29786 -- 2018-07-01 20:56:19
# 1214 of 24904: Fis-30623 -- 2018-07-01 21:15:14
# 24904 of 24904: WBD-Eup-82 -- 2018-07-02 05:36:41
i_all <- 24904; t_beg <- as_datetime("2018-07-01 20:49:01")
i_now <- 1214; t_now <- as_datetime("2018-07-02 05:36:41")
difftime(t_now, t_beg, units="hours") / i_now * (i_all - i_now) + t_now

# calculate avg file size ----
dir_spp <- "/Volumes/Best HD/mbon_data_big/aquamaps/species_rasters"
file_info(dir_ls(dir_spp)) %>% 
  summarize(
    size = mean(size)) %>%
  pull(size) * 24904 / 1000000

Load into PostGIS DB

gdal2wktraster.py -r C:\Temp\TutData\SRTM\tif\*.tif -t srtm_tiled 
      -s 4326 -k 60x60 -I > C:\Temp\TutData\SRTM\srtm.sql

Range Normalized Cumulative Species Map

To demonstrate the methodology of “range normalization”, let’s compare a cumulative species map with and without range normalization for two species with very different ranges:

• contracted: North Atlantic right whale (Eubalaena glacialis)

• expansive: blue whale (Balaenoptera musculus)

For a given cell (\(c\)), we might simply take the sum of AquaMap’s Relative Environmental Suitability (RES; \(e\); ranging in value [0, 1]) per species (\(s\)) across the number of species (\(S\)) to arrive at a cumulative species richness (\(R\)) that is weighted by the species suitability within the given cell:

\[ {R_c} = \sum_{s=1}^{S} e_{s,c} \]

This approach however gives equal weight for any given cell to a species with a contracted range as an expansive one. From a population perspective, this overrepresents the expansive range. To account for small areas with high degrees of unique endemism, we can “range normalize” by dividing the given relative environmental suitability for the species and cell (\(e_{s,i}\)) by the sum of the species relative environmental suitability (\(E_s\)) across all cells (\(C\)), i.e. the whole range, and further account for cell area (\(a_c\)) and IUCN extinction risk as a species weight (\(w_s\)):

\[ E_s = \sum_{c=1}^{C} e_{s,c} * a_c \\ {RN\ R_c} = \sum_{s=1}^{S} \frac{ e_{s,c} * a_c * w_s }{ E_s } \]

• \(a_c\): area of cell
• \(c\): cell
• \(e\): AquaMap’s Relative Environmental Suitability (RES; ranging in value [0, 1])
• \(E_s\): total environmental suitability, per species
• \(R_c\): species richness, per cell
• \(RN R_c\): range normalized species richness, per cell
• \(s\): species
• \(S\): total number of species
• \(w_s\): IUCN extinction risk weight of species

Note that the (IUCN extinction risk) weight (\(w_s\)) is not part of the range normalization denominator (\(E_s\)) so when summing species those with higher extinction risk are given more weight.

# species id's for right whale and blue whale
sid_blue  <- "ITS-Mam-180528"
sid_right <- "ITS-Mam-180537"

# IUCN extinction risk status
ext_blue  <- "EN"
ext_right <- "CR"

# extinction risk weights
# bumping up LC to 0.2 to give credit for presence and not expecting extinct 1.2
ext_wts <- c(LC=0.2, NT=0.4, VU=0.6, EN=0.8, CR=1, EX=1.2)
w_blue  <- ext_wts[ext_blue]
w_right <- ext_wts[ext_right]

get_am_raster <- function(sid, w=NULL, weight_area=T, range_normalize=T, attr_leaflet=T, extend_global=T){
  # Get AquaMaps raster
  # TODO: threshold param for binary output
  # sid: species id # sid <- sid_right
  # sid=sid_right; weight_area=T; range_normalize=T; attr_leaflet=T; extend_global=T; w=w_right
  
  # read tif
  tif <- glue("{dir_spp}/{sid}.tif")
  r <- raster(tif, crs=leaflet:::epsg4326)
  
  # set boundaries
  b <- extent(r)
  
  a = "attributes" # names(attributes(r)) # names(attributes(a))
  attr(a, "bounds") <- list(lng1 = b@xmin, lat1 = b@ymin, lng2 = b@xmax, lat2 = b@ymax) 

  # extend to global for calculations
  r <- extend(r, r_na)
  
  attr(a, "cell_avg_e") = cellStats(r, "mean", na.rm=T)
  attr(a, "sum_e") <- cellStats(r, stat='sum', na.rm=T)
  
  if (weight_area){
    r <- r * r_a
    
    attr(a, "cell_avg_e_a") = cellStats(r, "mean", na.rm=T)
    attr(a, "sum_area_km2")   <- cellStats(mask(r_a, r), stat='sum', na.rm=T)
    attr(a, "sum_e_area_km2") <- cellStats(r, stat='sum', na.rm=T)
  }
  
  if (range_normalize){
    E_s <- ifelse(weight_area, attr(a, "sum_e_area_km2"), attr(a, "sum_e"))
    r <- r / E_s
    
    attr(a, "E_s") <- ifelse(weight_area, "sum_e_area_km2", "sum_e")
    attr(a, "cell_avg_e_a_Es") = cellStats(r, "mean", na.rm=T)
  }
  
  if (!is.null(w)){
    r <- r * w
    
    attr(a, "w") <- w
    attr(a, "cell_avg_w") = cellStats(r, "mean", na.rm=T)
  }
  
  attr(a, "cell_avg") = cellStats(r, "mean", na.rm=T)
  
  if (attr_leaflet){
    r_leaflet <- projectRasterForLeaflet(raster_trim(r), method="ngb")
    attr(a, "r_leaflet") <- r_leaflet 
  }
  
  if (!extend_global)
    r <- raster_trim(r)

  # append attributes  
  walk(names(attributes(a)), function(x){
    attr(r, x) <- attributes(a)[[x]]
    assign("r", r, pos = parent.frame(n=2), inherits = FALSE)})
  #names(attributes(a))
  #names(attributes(r))
  
  r
}

# rasters for species relative environmental suitability [0 - 1]
r_right <-  get_am_raster(sid_right, weight_area=F, range_normalize=F) # plot(r_right)
r_blue  <-  get_am_raster(sid_blue, weight_area=F, range_normalize=F)   # plot(r_blue)
r_right_w_a_r <-  get_am_raster(sid_right, w=w_right, weight_area=T, range_normalize=T) # plot(r_right)
r_blue_w_a_r  <-  get_am_raster(sid_blue, w=w_blue, weight_area=T, range_normalize=T)   # plot(r_blue)

# names(attributes(r_right))
# attr(r_blue, "cell_avg")
# attr(r_right, "cell_avg")
# attr(r_right, "sum_area_km2")

sum_am_rasters <- function(..., attr_leaflet=T){
#add_am_rasters <- function(..., spp_wts=NULL, spp_range_normalize=T, spp_all_normalize=T){
  # spp_range_normalize: for each species map, divide by sum of all pixels, before summing
  # spp_all_normalize: after summing all species maps, divide by max pixel value to get range [0,1]
  # spp_wts: species weights, to apply in same order as rasters. defaults to 1
  # TODO: add species richness factor, so more species -> stronger result
  # ...: rasters
  
  rasters <- list(...)
  s <- stack(rasters)
  
  # if (!is.null(spp_wts)){
  #   stopifnot(length(rasters)==length(spp_wts))
  #   for (i in 1:length(spp_wts))
  #     s[[i]] <- raster(s, layer=i) *  spp_wts[i]
  # } #else {
    #spp_wts = rep(1, length(rasters))
  #}
    
  # spp_wts <- ext_wts[c(ext_right, ext_blue)]
  
  # get weighted mean
  #r <- weighted.mean(s, spp_wts)
  
  # sum across all species rasters
  #r <- calc(s, function(x) )
  #r <- calc(s, fun = sum, na.rm = T)
  r <- sum(s, na.rm = T)
  #plot(raster(s, layer=1), col=rev(rainbow(255)))
  
  # mask
  m <- overlay(s, fun=function(x) ifelse(sum(!is.na(x)) > 0, 1, NA))
  r <- mask(r, m)
  
  # trim
  #r <- 
  #summary(r); plot(r, col=rev(rainbow(255)))
  #summary(r0); plot(r0, col=rev(rainbow(255)))
  #summary(r0); plot(stack, r0)
  #r <- calc(s, fun = sum, na.rm = T)
  #summary(r0)
  
  # if (spp_all_normalize)
  #   r <- r / cellStats(r, stat='max', na.rm=T)
  
  # project for leaflet
  if (attr_leaflet)
    attr(r, "r_leaflet")  <- projectRasterForLeaflet(raster_trim(r), method="ngb")
  
  r
}

#r_sum   <- add_am_rasters(r_right, r_blue, spp_range_normalize=F, spp_all_normalize=F, spp_wts = NULL)
#r_nrsum <- add_am_rasters(r_right, r_blue, spp_range_normalize=T, spp_all_normalize=F, spp_wts = NULL)
r_sum   <- sum_am_rasters(r_right, r_blue) # plot(r_sum)
r_nrsum <- sum_am_rasters(r_right_w_a_r, r_blue_w_a_r)

#plot(r_sum)
#plot(r_nrsum)
#r_add3 <- add_am_rasters(r_right, r_blue, spp_range_normalize=T, spp_all_normalize=T, spp_wts = NULL)
#r_add4 <- add_am_rasters(r_right, r_blue, spp_range_normalize=T, spp_all_normalize=T, 
#                         spp_wts = ext_wts[c(ext_right, ext_blue)])

# summary(r_sum)
# summary(r_nrsum)
# summary(r_add3)
# summary(r_add4)

fit_am_raster <- function(map, r){
  bounds <- attr(r, "bounds")
  with(bounds, fitBounds(map, lng1, lat1, lng2, lat2))
}

add_am_raster <- function(m, r, lyr_name, r_type="spp"){
  # r_type="spp" or "sum"
  # r <- r_right
  r <- attr(r, "r_leaflet")
  
  max_r <- cellStats(r, stat='max', na.rm=T)
  sum_r <- cellStats(r, stat='sum', na.rm=T)

  rng <- c(0, 1)
  if (r_type != "spp"){
    rng <- c(0, max_r)
  }
  
  pal <- colorNumeric(
    "inferno", rng, na.color = "transparent", reverse=T)

  lab_sci <- function(v, type="numeric"){
    #browser()
    if (r_type == "spp"){
      labs <- case_when(
        v == 0 ~ "0, 0",
        v == last(v) ~ as.character(glue(
          "{v}, {formatC(v/sum_r, format = 'e', digits = 1)}")),
        TRUE ~ "")
    } else {
      #if (lyr_name == "sum NR RES") browser()
      labs <- case_when(
        v == 0 ~ "0, 0",
        v == last(v) ~ as.character(glue(
          "{v}, 1")),
        TRUE ~ "")
    }
    labs
  }
  
  m %>%
    addRasterImage(
      r, colors = pal, opacity = 0.6, project=F, 
      group=lyr_name) %>%
    addLegend(
      group = lyr_name, pal = pal,
      position = "bottomright",
      values = rng, labFormat = lab_sci,
      title = ifelse(
        r_type == "spp",
        glue("{lyr_name} (RES, E_s)<br>&nbsp;&nbsp;sum={format(sum_r, big.mark=',', digits=1, nsmall=1)}"),
        glue("{lyr_name} (max, 1)")))
}

map_am_layer <- function(lyr){
  
  m <- leaflet() %>%
    addProviderTiles(
      providers$Esri.OceanBasemap, options=providerTileOptions(opacity=0.6)) %>%
    add_am_raster(r_right, "right whale") %>%
    add_am_raster(r_blue, "blue whale") %>%
    add_am_raster(r_sum, "sum", r_type="sum") %>%
    add_am_raster(r_nrsum, "sum RN", r_type="sum") %>%
    addLayersControl(
      baseGroups = c("ESRI Ocean Basemap"),
      overlayGroups = c("right whale", "blue whale", "sum", "sum RN"),
      #overlayGroups = c("right whale", "blue whale", "sum"),
      options = layersControlOptions(collapsed = T)) %>%
    fit_am_raster(r_blue) %>%
    addMiniMap(position = "bottomleft", toggleDisplay = T, minimized = T) %>%
    hideGroup("blue whale") %>%
    hideGroup("right whale") %>%
    hideGroup("sum") %>%
    hideGroup("sum RN") %>%
    showGroup(lyr)
  m
}

map_am_layer("blue whale")

map_am_layer("right whale")

map_am_layer("sum")

map_am_layer("sum RN")

Range Normalization, Regardless of Cell Area

The equation above is problematic because it biases equatorial cells and downweights poleward cells. Another approach is to ask what does “1” mean per species richness factor?

• Extinction risk (\(w_s\)): 1 = extinct, 0.8 = critically endangered, …, 0.2 = least concern

• Range normalization (\(N_s\)): species normalization factor, based on the smallest species range within a given taxonomic group, eg Mammalia.

So alternatively, we describe a range normalization factor per species (\(N_s\)) that can be applied to all cells (\(C\)) regardless of the cell area (ie largest at equator, smallest at pole) relative to some minimum range size we define as 1.

\[ E_s = \sum_{c=1}^{C} e_{s,c} * a_c \\ N_s = \frac{min\{E_1, ..., E_S\}}{E_s} \\ NR_c = \sum_{s=1}^{S} e_{s,c} * N_{s} * w_s \]

• \(E_s\): sum of environmental suitability across range, ie all cells, per species \(_s\)
• \(e_{s,c}\): AquaMap’s Relative Environmental Suitability (RES; ranging in value [0, 1]) per species \(_s\) and cell \(_c\)
• \(a_c\): area of cell \(_c\)
• \(C\): all cells
• \(N_s\): range normalization factor, per species \(_s\)
• \(\{E_1, ..., E_S\}\): set of species environmental suitability sums \(E_s\) across all species \(_S\)
• \(NR_c\): range normalized species richness per cell \(_c\), also weighted by extinction risk \(w\)
• \(w_s\): IUCN extinction risk weight of species

New Range Normalized Mammalia Richness

Get species in Class “Mammalia” and IUCN extinction risk…

# species id's for right whale and blue whale
# sid_blue  <- "ITS-Mam-180528"
# sid_right <- "ITS-Mam-180537"

library(scales)
library(rredlist)
select = dplyr::select

get_iucn_category <- function(genus_species){
  # token issued to ben@ecoquants.com via http://apiv3.iucnredlist.org/api/v3/token
  rl_token <- "fb963555580923eb0f0b2778060e5bcbd0cca50c0b21c416ed7ff68da9c25c11"
  #genus_species <- spp_g$genus_species[1]
  
  #genus_sp <- 'Balaenoptera brydei'
  #genus_sp <- "Mesoplodon hectori"
  #genus_sp <- "Mesoplodon hectori"
  
  r <- rl_search(genus_species, key=rl_token)
  # r <- rl_search_(name=genus_species, id=NULL, region=NULL, key=rl_token)
  # r <- rredlist:::rr_GET(rredlist:::nir("species", "species/id", name=genus_species, id=NULL, region=NULL), id=NULL, region=NULL, key=rl_token)
  # r <- rredlist:::nir("species", "species/id", name=genus_species, id=NULL, region=NULL)
  #browser()
  if (is_empty(r$result)){
    iucn_cat <- NA
  } else {
    iucn_cat <- r$result$category
  }
  #cat(glue("{genus_species}: {iucn_cat}"), "\n")
  
  iucn_cat
}

# bumping up LC to 0.2 to give credit for presence and not expecting extinct 1.2
iucn_weights <- c(LC=0.2, NT=0.4, VU=0.6, EN=0.8, CR=1, EX=1.2)

redo_spp_g <- F
if (!file.exists(spp_csv) | redo_spp_g){
  spp_g <- spp %>%
    filter(Class == "Mammalia") %>%
    arrange(Genus, Species) %>% 
    mutate(
      group         = "Mammalia",
      genus_species = glue("{Genus} {Species}") %>% as.character(),
      iucn_category = map_chr(genus_species, get_iucn_category),
      iucn_weight   = iucn_weights[iucn_category],
      tif           = glue("{dir_spp}/{SPECIESID}.tif"),
      stack_name    = str_replace_all(SPECIESID, "-", "."))
  
  write_csv(spp_g, spp_csv)
}
spp_g <- read_csv(spp_csv)

# get raster stack of species, extended to globe
s_g <- stack(sapply(spp_g$tif, function(tif) raster(tif, crs=leaflet:::epsg4326) %>% extend(r_na)))
names(s_g) <- spp_g$SPECIESID

get_spp_stat <- function(stack_name, s, stat="sum_RES"){
  #stack_name <- spp_g$r[1] 
  r <- raster(s, stack_name)
  if (stat == "sum_RES"){
    cellStats(r, stat = "sum")
  } else if (stat == "sum_area_km2"){
    return(cellStats(mask(r_a, r), stat = "sum"))
  }
}

# get stats on sum RES and area
spp_g <- spp_g %>%
  mutate(
    sum_RES        = map_dbl(stack_name, get_spp_stat, s_g, "sum_RES"),
    sum_area_km2   = map_dbl(stack_name, get_spp_stat, s_g, "sum_area_km2"),
    rescale_RES    = rescale(1/sum_RES),
    logrescale_RES = rescale(1/log(sum_RES)))

write_csv(spp_g, spp_csv)
spp_g <- read_csv(spp_csv)

spp_g %>%
  select(SPECIESID, FBname, genus_species, iucn_category, iucn_weight, sum_area_km2, sum_RES, rescale_RES, logrescale_RES) %>%
  DT::datatable()

table(spp_g$iucn_category, useNA = "ifany")

## 
##   CR   DD   EN   EX   LC   NT   VU <NA> 
##    2   39   13    1   47    4   10    3

iucn_weights

##  LC  NT  VU  EN  CR  EX 
## 0.2 0.4 0.6 0.8 1.0 1.2

TODO: Range normalize across Mammalia based on latest equation using minimum species range as reference point.

library(DT)

hist(spp_g$sum_RES)

hist(spp_g$rescale_RES)

hist(spp_g$logrescale_RES)

with(spp_g, plot(sum_RES, rescale_RES))

with(spp_g, plot(sum_RES, logrescale_RES))

sum_am_stack <- function(s, spp, do_rel=T, do_ext=T, do_rng=T, do_sca=T, binary_threshold = 0.4, log_rng=T){

  # parameters:
  # s = # stack of species AM RES rasters for group
  # spp = # table of species information
  # do_rel = T # relative environmental suitability; otherwise binary
  # binary_threshold = 0.4 # if do_ext=F
  # do_ext = T # extinction risk
  # do_rng = T # range contraction, ie endemism
  # do_sca = T # final rescale [0-1]
  
  # do_res: relative environmental suitability, otherwise binary
  if (!do_rel){
    v_to_binary <- function(v, na.rm=T){
      if (is.na(v)) return(NA)
      ifelse(v >= binary_threshold, 1, NA)
    }
    
    s <- stackApply(s, 1:nlayers(s), v_to_binary)
    names(s) <- names(s)
  }
  
  # do_ext: extinction risk
  if (do_ext){
    s <- s * spp$iucn_weight
  }
  
  # do_rng: range contraction, ie endemism
  if (do_rng){
    if (log_rng){
      s <- s * spp$logrescale_RES
    } else {
      s <- s * spp$rescale_RES
    }
  }
  
  # sum
  r <- sum(s, na.rm = T)
  
  # do_sca: final rescale [0-1]
  if (do_sca){
    values(r) <- rescale(values(r))
  }
  
  # mask
  m <- overlay(s, fun=function(x) ifelse(sum(!is.na(x)) > 0, 1, NA))
  r <- mask(r, m) # plot(r); 
  # projectRasterForLeaflet(raster_trim(r), method="ngb") %>% mapview::mapview(attr(r, "r_leaflet"))

  r
}

r_bes  <- sum_am_stack(
  s_g, spp_g, do_rel=F, do_ext=T, do_rng=F, do_sca=T, binary_threshold = 0.4)
r_be  <- sum_am_stack(
  s_g, spp_g, do_rel=F, do_ext=T, do_rng=F, do_sca=F, binary_threshold = 0.4)
r_re  <- sum_am_stack(
  s_g, spp_g, do_rel=T, do_ext=T, do_rng=F, do_sca=F)
r_res  <- sum_am_stack(
  s_g, spp_g, do_rel=T, do_ext=T, do_rng=F, do_sca=T)
r_regs <- sum_am_stack(
  s_g, spp_g, do_rel=T, do_ext=T, do_rng=T, log_rng=F, do_sca=T)
r_rels <- sum_am_stack(
  s_g, spp_g, do_rel=T, do_ext=T, do_rng=T, log_rng=T, do_sca=T)

map_am_sum <- function(r, name){
  r_m <- projectRasterForLeaflet(raster_trim(r), method="ngb")
  
  rng <- range(values(r), na.rm=T)
  pal <- colorNumeric(
    "inferno", rng, na.color = "transparent", reverse=T)

  leaflet() %>%
    addProviderTiles(providers$Stamen.TonerLite) %>%
    addRasterImage(
      r_m, colors = pal, opacity = 0.8, project = F, group = name) %>%
    addLegend(
      group = name, pal = pal,
      position = "bottomright",
      values = rng, 
      title = name)
}
map_am_sum(r_bes,  "Mammalia: Binary<br> * Extinction")

map_am_sum(r_res,  "Mammalia: Relative<br> * Extinction")

map_am_sum(r_regs, "Mammalia: Relative<br> * Extinction<br> * Endemism")

map_am_sum(r_rels, "Mammalia: Relative<br> * Extinction<br> * log(Endemism)")

Observation Coverage

Compare species ranges with available OBIS data to determine degree of coverage, especially for determining where we expect concentrations of endangered species and yet have no observational data.

For example with North Atlantic Right Whale (Eubalaena glacialis):

• OBIS observations: 5695 records for Eubalaena glacialis

  

• AquaMaps Relative Env Suitability

• IUCN Range Map

  

The topic of species home ranges is rich. Am keeping track of references in the mbon-indicators Zotero group, eg folder species home ranges:

• Home range - Wikipedia

• adehabitatHR.pdf

• Burgman, M. A., & Fox, J. C. (2003). Bias in species range estimates from minimum convex polygons: implications for conservation and options for improved planning. Animal Conservation, 6(1), 19–28. https://doi.org/10.1017/S1367943003003044 Calenge, C. (n.d.). Home Range Estimation in R: the adehabitatHR Package, 61.
• Elith, J., & Leathwick, J. R. (2009). Species Distribution Models: Ecological Explanation and Prediction Across Space and Time. Annual Review of Ecology, Evolution, and Systematics, 40(1), 677–697. https://doi.org/10.1146/annurev.ecolsys.110308.120159
• Gaston, K. J., & Fuller, R. A. (2009). The sizes of species’ geographic ranges. Journal of Applied Ecology, 46(1), 1–9. https://doi.org/10.1111/j.1365-2664.2008.01596.x
• Getz, W. M., Fortmann-Roe, S., Cross, P. C., Lyons, A. J., Ryan, S. J., & Wilmers, C. C. (2007). LoCoH: Nonparameteric Kernel Methods for Constructing Home Ranges and Utilization Distributions. PLoS ONE, 2(2), e207. https://doi.org/10.1371/journal.pone.0000207
• Getz, W. M., & Wilmers, C. C. (2004). A local nearest-neighbor convex-hull construction of home ranges and utilization distributions. Ecography, 27(4), 489–505. https://doi.org/10.1111/j.0906-7590.2004.03835.x
• Graham, C. H., & Hijmans, R. J. (2006). A comparison of methods for mapping species ranges and species richness. Global Ecology and Biogeography, 15(6), 578–587. https://doi.org/10.1111/j.1466-8238.2006.00257.x
• Kearney, M., & Porter, W. (2009). Mechanistic niche modelling: combining physiological and spatial data to predict species’ ranges. Ecology Letters, 12(4), 334–350. https://doi.org/10.1111/j.1461-0248.2008.01277.x

• Peterson, A. T. (2017). Problems with reductive, polygon-based methods for estimating species’ ranges: reply to Pimm et al. 2017: Estimating Species’ Ranges. Conservation Biology, 31(4), 948–951. https://doi.org/10.1111/cobi.12929

Research